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Abstract 



■ The zero energy modes induced by vacancies in ABC stacked trilayer graphene are investigated. Depending on the position of 
^ ]the vacancy, a new zero energy solution is realised, different from those obtained in multilayer compounds with Bernal stacking. 
&DThe electronic modification induced in the sample by the new vacancy states is characterised by computing the local density of 
^ "states and their localisation properties are studied by the inverse participation ratio. We also analyse the situation in the presence of 
a gap in the spectrum due to a perpendicular electric field. 
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1. Introduction 

Recent experimental advances aiming to generate better sam- 
ples for electronic devices have allowed the obtention of high 
quality samples not only of monolayer graphene but also of bi- 
layer graphene (BLG) and trilayer graphene (TLG) man. 
One of the major problems preventing applications of mono- 
layer graphene is the difficulty to open and control a gap in 
the samples. To this respect bilayer and multilayer samples are 
more promising J]]]. The band structure of few layer graphenes 
depends on the stacking order |4l|5j,|6|], what offers the possibil- 
ity of tuning electronic properties. After the great excitement 
awaken by the bilayer compound due to fhepossibility of open- 
ing a tunable gap with an external gate 10,18, 01 me interest has 
moved recently to the trilayer materials due to their enigmatic 
properties. As in the Bernal AB stacked BLG an external elec- 
tric field allows a tunable band gap in the ABC-stacked TLG 
"S HI HI while the ABA-stacked TLGs are semimetals with 
electric field tunable band overlap |2, 121. 

Structural defects - vacancies, ad-atoms, and other - which 
may appear during the fabrication process are very important in 
the graphene materials. While in the first times after the synthe- 
sis the concern was that they may deteriorate the performance 
of graphene-based devices, later tendencies point to their pos- 
itive use in some applications, as they make it possible to tai- 
lor the local properties of graphene and to achieve new func- 
tionalities [13]. Being one-atom thick the graphene materials 
are extremely sensitive to the presence of adsorbed atoms and 
molecules and, more generally, to defects such as vacancies, 
holes and/or substitutional dopants. This property, apart from 
being directly usable in molecular sensor devices, can also be 
employed to tune graphene electronic properties. The possi- 
bility of a controlled manipulation of atoms and molecules on 
graphene has opened a new area of research that allows to ob- 
serve chemical interactions or structural modifications of low 
contrast molecules or nano-objects J14I1 . Among the defect- 
induced properties in graphene, magnetism is one of the most 



appealing 1 15, Hit 17 , lUl . Recent experimental findings [19] 
have ascertain the role of vacancies in the magnetic properties 
of the material. 

Vacancies are recognised as important scattering centers in 
monolayer and BLG 12011 and are known to induce zero-energy 
modes 12111 that modify the low-energy properties of the sam- 
ples. The existence and nature of localised states arising from 
vacancies in BLG was analyzed in a recent paper [22]. It 
was found that the two different types of vacancies that can be 
present in the bilayer system give rise to two different types of 
states: quasilocalised states, decaying as 1/r at long distances r 
from the vacancy, similar to these found in the monolayer case 
[23], and truly delocalised states whose wave function remains 
constant as r — > 00. These later new midgap states were found 
to become localised inside the gap induced in the bilayer sam- 
ple by the electric field effect. The analysis of these vacancy- 
induced states was generalized in [22] to multilayer graphene 
systems with ABAB- • • Bernal stacking. The recent experi- 
mental and theoretical activity around the trilayer compounds 
il|I|25l|2l|27|] has shown that the system in the rhombohe- 
dral ABC stacking has very different electronic properties than 
its Bernal partner. 

In this paper we will revise the situation of the midgap 
states induced by vacancies in few layers graphene, with espe- 
cial attention to TLG. We will see that the ABC-rhombohedral 
stacked presents yet a new type of zero-energy state with no 
analogue in the mono or bilayer compounds. The characterisa- 
tion of this new state exhausts the possibilities for the vacancy 
states in multilayer graphene with the usual stacking. The pa- 
per is organized as follows: in section|2]we give a summary of 
the situation encountered in the bilayer system and its exten- 
sion to multilayer graphene with Bernal (AB) stacking. Next 
we explain in section [3] the new aspects encountered in the tri- 
layer material with rhombohedral (ABC) stacking. In section 
[4] we summarise the analysis of the work and discuss possible 
physical consequences. 
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Figure 1: (color online), (a) Bilayer lattice structure and main tight-binding pa- 
rameters, (b) The general shape of the density of states for the minimal model. 
The inset shows the changes in the DOS near the Fermi point when the various 
tight binding parameters are included. 



2. Vacancy states in Bernal stacked multilayer graphene 

The study and characterisation of vacancy states in mono- 
layer graphene has been an important subject that started prior 
to the synthesis of the material and continues to our days 
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driven in part by the search for magnetism in sp car- 
bon compounds 111 611 . In the monolayer case there is only one 
type of vacancy whose presence induces two degenerate quasi- 
localised states at the Fermi energy. An analytic construction of 
the wave function of the vacancy was done in 112.311 by matching 
surface state solutions at zigzag edges with those localised at 
Klein edges for a suitable boundary condition. In the contin- 
uum limit the wave function can be written as 



,/K'.r 



x + ly x — ly 



(1) 



where K and K' are the reciprocal space vectors of the two 
inequivalent corners of the first Brillouin zone, and (x,y) are 
distances in a reference frame centred at the vacancy position. 
The wave function is peaked at the position of the vacancy and 
decays as 1 jr away from it, a behaviour termed quasi-localised. 

This procedure to describe vacancy states was generalised to 
the AB stacked BLG in [22]. Because of the similitudes with 
the trilayer case and to fix the notation we will present the main 
features of the bilayer in some detail in what follows. 

The lattice structure of a BLG is shown in Fig. Q] In the 
AB-Bernal stacking the top layer is shifted with respect to the 
bottom layer by one C-C distance. As a result only half of the 
atoms in any layer have a direct neighbour joined by j\ in the 
other layer. We use indices 1 and 2 to label the top and bottom 
layer, respectively. The main tight-binding hopping parameters 
are also shown in the figure as well as the different electronic 
structure near the Fermi level for different values of the tight 
binding parameters. In the minimal model adopted in most of 
the works only j\ is different from zero [7]. The estimated 
values of the parameters are t ~ 3eV and y\ jt ~ 0.1 lfl5ll . 

In the Bernal BLG there are two different types of vacancies 
giving rise to unpaired atoms: vacancies located at A1/B2 or 
B1/A2 [see Fig. [TJa)]. The first type A1/B2 is produced by 
removing a site having a neighbour in the adjacent layer and 
is usually named B vacancy. The second type BljAl resulting 
when the removed site is not connected to the other layer is 
called a. 

We obtained an analytical solution for the states induced by 
these two types of vacancies in 12211 generalising the procedure 



done in the monolayer case [ 23] . For vacancies at A1/B2 sites 
the wave function obtained is the same given by Eq. ([]]) which 
corresponds to a quasilocalised, zero-energy state decaying as 
1 jr around the vacancy in the same layer where the vacancy sits 
but in the opposite sublattice. This vacancy has then the same 
properties as these found in the monolayer case. The new state 
found in 12211 corresponds to vacancies at B1/A2 lattice sites. 
The solution in this case has the form 



r(x,y) 



t x + iy 



t x — ly 



, (2) 



where T'fx, y) is the quasi-localised state given in Eq. ([TJ, and 
the two component wave function T ~ [<pu <p2\ refers to the 
two layers; first and second components for the first and second 
layers, respectively. This is a delocalised state, with the pecu- 
liarity of being quasi-localised in one layer (where the vacancy 
sits) and delocalised in the other where it goes to a constant 
when r — > oo. 

The same analytical construction followed for the minimal 
model in BLG can be directly applied to multilayer graphene in 
Bernal stacking and for the particular case of the ABA TLG 
The quasi-localised state ((TJ is a solution in any multilayer 
with a Al/B2-vacancy. For a fil/A2-vacancy the solution is 
a generalisation of state (f2| with a quasi-localised component 
in the layer where the vacancy resides and delocalised compo- 
nents in the layers right on top and below this one: <b(x,y) ~ 
[^¥(x,y), T2(x,y), Y2(x,y)], where T2(x,y) refers to the second 
component in the right hand side of Eq. (0. 

The behaviour of the vacancy states was ascertained in the 
bilayer system by numerical calculations. A summary of the 
findings of our previous work on Bernal stacked multilayer 
graphene is the following: Associated to the presence of va- 
cancies and to the existence of a gap in the spectrum gener- 
ated by an electric field E z three different types of behaviour for 
vacancy-induced states are found: 

1. For E z = a B— vacancy (A1/B2 site) induces quasi- 
localized state (1/r behaviour). 

2. For E z = an a-vacancy (B 1 /A2 site) induces a resonance 
due to a delocalized state. 

3. For E z + a B— vacancy produces a resonance inside the 
continuum near the band edge while an a-vacancy gives 
rise to a truly localized state inside the gap. This is the 
most interesting state for the magnetic implications. 



3. ABC trilayer 

3.1. Model 

We follow the tight-binding approximation and consider the 
minimal model where the in-plane hopping energy, t, and the 
inter-layer hopping energy, j\, define the most relevant energy 
scales. This is the minimal model for AB stacked multilayer 
graphene described above. The simplest tight-binding Hamilto- 
nian describing non-interacting 7r-electrons in TLG then reads: 
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Figure 2: (color online) (a) Trilayer lattice structure, (b) Scheme indicating 
in-plane (full line) and inter-layer (dashed line) nearest neighbor hopping be- 
tween different sublattices. (c) LDOS at the three non-equivalent sites in perfect 
trilayer. 



Vacancies are modelled as missing sites in Eq. (0. There are 
three non-equivalent vacancy sites. As can be inferred from the 
lattice sketch shown in Fig.[2jb), these sites correspond to a va- 
cancy occurring in sublattice A1/B3, A2/B2, and A3/B1. The 
LDOS at these three lattice sites is shown in Fig. |2jc) for the 
perfect lattice (no vacancy present). Despite the oscillatory be- 
haviour due to finite size effects, distinct thermodynamic limit 
features can be seen |,29?1. The low energy physics is determined 
by sublattices A3/B1, and at E = a Van-Hove singularity de- 
velops due to the cubic spectrum E ~ k 3 . Sublattices A1/B3 
and A2/B2 have a vanishing contribution to the density of states 
at low energy. Due to interlayer hybridization their contribution 
is appreciable only for \E\ > j\. 



1=1 R,cr 

(3) 

with H, being the SLG Hamiltonian 

Hi = -t [a V(R) + aljR)bi,AR - a, ) + (4) 
a|; cr (R)V(R-a 2 ) + h.c] > 

where a !j0 -(R) [£>, j0 -(R)] is the annihilation operator for electrons 
at position R in sublattice Ai (Bi), i = 1, 2, 3, and spin <x, and 
ai and a2 are the primitive vectors of the underlaying Bravais 
lattice. If a perpendicular electric field is applied 111 ill , the fol- 
lowing potential energy should be added to Eq. (0, 



(5) 



with n, jCr = a". (R)a, ](r (R) + b". (R)Z> !j0 -(R). We use the same 
values for in-plane and inter-layer hopping as for AB stacked 
multilayer graphene IU5I1 . which, as we have seen above, im- 
ply yi/t ~ 0.1 «C 1. Extra hopping terms, as long as they 
preserve the bipartite nature of the lattice, may introduce quan- 
titative changes but not qualitative 11231 l22l 12811 . Hopping terms 
that break the bipartite character of the lattice may be treated 
perturbatively afterwards [18]. 

The main feature of the ABC compound that makes it dif- 
ferent from its Bernal ABA counterpart is the lack of mirror 
symmetry with respect to the middle layer. As discussed in 
Ir29l l30ll this type of staking allows the derivation of a low 
energy effective hamiltonian which involves only the unlinked 
atoms (Ai, B3) given by 



qjeff = _ V J_ 
2 

n 





k 3 



k* 3 




(6) 



where k = k x + ik v and Vf = 3ta/2, with a the C-C distance. It 
also guarantees the topological stability of the low-energy chiral 
effective Hamiltonian and makes its behavior similar to the AB 
bilayer discussed previously [30]. As it happens in the bilayer 
case, a gap can be induced in the ABC system by an external 
gate. 



3.2. New midgap state: Analytical considerations 

From what is known for the single layer and for bernal 
stacked multilayer graphene we may immediately infer the fol- 
lowing. 

1 . For a vacancy at any of the two possible places in the mid- 
dle layer of the ABC TLG the zero energy mode is of the 
new bilayer type: quasi-localised in the middle layer, and 
delocalised in the adjacent layer not connected to the va- 
cancy. 

2. For a vacancy at any of the outer layers, and residing in 
the sublattice which is connected to the middle layer, then 
we simply have the single layer solution: a state quasi- 
localised around the vacancy. 

3. For a vacancy at any of the outer layers, but residing in the 
sublattice which is not connected to the middle layer we 
have a new solution with amplitude over the three layers 
simultaneously. 

By extension of the bilayer case we expect the new solution to 
have a quasi-localised component in the layer with the vacancy 
and a delocalised one in the middle layer. This is nothing but 
the new bilayer solution for these two components. The cor- 
rect behaviour of the third component may be determined by 
the matching procedure used previously lEM E2I1 . We do not 
follow such prescription here. Instead, we perform a numerical 
analysis detailed in the next section. Nevertheless, we may ex- 
pect this new component to be delocalised as well just based on 
the continuum, low energy model, Eq. (0 that describes ABC 
TLG. Far from the vacancy position such component must sat- 
isfy dltyiz, z) — 0. As a natural extension of the bilayer case we 
have the solution ifr(z,z) = z 2 f(z), with /(I) meromorphic. As- 
suming, as in the single layer and bilayer cases, that the mero- 
morphic function is the usual f(z) = 1 /z, the new component is 
is manifestly delocalised. 

3.3. Numerical analysis 

As mentioned before, vacancies are introduced in the TLG 
lattice by elimination of atom sites. In our approach we use 
a simple tight binding Hamiltonian and do not include any re- 
construction in the remaining structure. We have analysed the 
changes induced in the LDOS for sites around each vacancy 
for the three different types described. The LDOS is com- 
puted by the recursive Green's function method in clusters with 
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Figure 3: (color online) Upper: LDOS for a vacancy at sublattice (a) A1/B3, 
(b) A2/B2, and (c) A3/B1. The LDOS is computed at a lattice site closest to 
the vacancy. Down: IPR for the zero-energy mode induced by the previous 
vacancies. In panels (a) and (b) linear guide lines are shown to illustrate the 
scaling. In panel (c) the line is a fit with linear and quadratic terms. 

W = 3 x 1200 2 , from which the thermodynamic limit can be 
inferred. 

The localisation character of vacancy-induced modes is stud- 
ied through finite-size-scaling of the inverse participation ratio 
(IPR). The later is defined as 

N 

P Y = J] |^ v (0| 4 (7) 
t 

for the eigenstate v, where tp v (i) is its amplitude at site i. We 
perform exact diagonalization on small clusters with N up to 
3 x 82 2 sites. The IPR for extended, quasi-localized, and truly 
localized states scales distinctively with N [31]. While for ex- 
tended states we have P v ~ N~ l , for quasi-localized states the 
1/r decay implies P v ~ \og(N)~ 2 (consequence of the definition 
of the IPR in terms of normalized eigenstates). For localized 
wavefunctions the significant contribution to P v comes from 
the sites in which they lie, and a size independent P v shows up. 

3.3.1. Gapless case 

The LDOS and corresponding IPR are shown in Fig.|3]for the 
three types of vacancy states discussed in this work. The LDOS 
shown in the upper part is computed at a lattice site closest to 
the vacancy. We can see a sharp resonant peak for the first type 
of vacancy in panel [3j a) that corresponds to the state common 
to the monolayer. The wave function has amplitude only in 
the layer of the vacancy and is quasi-localised in the opposite 
sublattice. The 1/r decay of the wave function is apparent from 
the IPR shown in the down part|3la). 

The middle panel |3b) depicts the results obtained for the 
vacancy A2/B2. This type of vacancy, also existing in the bi- 
layer AB, is quasi-localised in one layer (where the vacancy 



Figure 4: (color online) Upper: LDOS for a vacancy at sublattice (a) A1/B3, 

(b) A2/B2, and (c) A3/B1 for a finite gap, V = O.lr. The LDOS is computed at 
a lattice site closest to the vacancy. Localized modes inside the gap are signaled 
by the arrow. Down: IPR for a vacancy at sublattice (a) A1/B3, (b) A2/B2, and 

(c) A3/B1 for a finite gap, V = O.lf. As usually done in regions of continuum 
DOS, the IPR in (a) is an average over states in an energy bin Ail = 0.3f around 
the gap-edge resonance shown in the upper part (a). In (b) and (c) the IPR is for 
the in-gap mode shown in the upper part. Lines are guides to the eyes. 

sits) and delocalised in the adjacent layer where it remains con- 
stant when r — > oo. This behaviour can be seen as a broader 
feature in the LDOS (upper part[3jb)) that can be attributed to 
the quasi-localized component in the layer where the vacancy 
sits (and thus the feature). This interpretation is fully corrobo- 
rated by the IPR scaling analysis shown in the down part[3jb). 

Panel [3]c) shows the behaviour of the new vacancy state 
found in this work. The absence of a distinct resonance around 
zero energy is a clear indication that the wave function of this 
vacancy has amplitude in all three layers. The feature associ- 
ated with the quasi-localised component in the layer where the 
vacancy resides is thus weaker than in the other cases, and is 
burried in the continuum background of band states. The results 
for the IPR (lower panel) [3ja) confirm the behaviour predicted 
by the general arguments done in Sec. 13.21 The IPR scales as 
N~ l as corresponds to an extended state, even though finite size 
effects are quite strong in this case (higher powers of N~ l are 
necessary to fit the data). 

3.3.2. Gapped case 

The ABC TLG presents, as the AB BLG, the possibility of 
opening and controlling a gap in the spectrum by applying an 
external electric field E z . We have studied the behaviour of 
the midgap states when a gap opens. We consider a gap in- 
duced through a perpendicular electric field E z = V/ (ed), where 
d ~ 0.34 nm is the interlayer distance. Its presence is modelled 
by adding an on-site energy term: -V/2 at layer 1 and V/2 at 
layer 3, as given by Eq. ©. 

The results for the gapped case are shown in Fig. 13.3.21 We 
plot the LDOS (upper panel) and IPR (lower panel) for the three 
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different types of vacancies. 

The first type of vacancy state 13.3.21 a). corresponds to the 
monolayer type having amplitude only in the layer where the 
vacancy sits. From the LDOS it can be seen that the quasi- 
localised state shown in the upper panel of Fig. [3j a), becomes 
a resonance around +V/2 in the gaped case. The "monolayer" 
vacancyOa) becomes delocalised at the gap edge. The IPR for 
this case is an average over the gap-edge resonance shown in 
the LDOS. 

The LDOS for the vacancies of the types [3fb) and[3]c) (up- 
per panel) show that they live inside the gap. We have ascer- 
tained this result performing calculations for different gap sizes. 
Their asymmetric weight over the two layers explains why they 
appear off zero-energy. The IPR scaling to a constant (lower 
panel) shows that these vacancies give rise to truly localised 
states inside the gap. This is the same behaviour encountered 
for the bilayer case in I22I1 . 

4. Conclusions and discussion 

We have analysed the various midgap states that arise from 
the presence of vacancies in multilayer graphene with special 
emphasis on their degree of localisation. Previous works based 
on a generalisation of the construction of vacancy states in bi- 
layer AB to multilayer compounds with Bernal (AB) staking 
found very different localisation properties for the two differ- 
ent types of vacancies that are present in these compounds. 
We have completed these works by analysing the rombohedral 
ABC trilayer graphene and found yet another new type of va- 
cancy states. Based in tight binding arguments ane can show 
that the new vacancy state survives beyond three layers. In the 
rhombohedral multilayer the amplitude spreads over all layers, 
though decaying exponentially as {y\ / ?)" away from the layer 
where the vacancy sits (n = 0). We can also apply continuum 
arguments. It was shown in 1 30] that the general rombohedral 
staking including the links (A\ - ^2,^2 - #3, ■ ■ ■ ,A/v-i - Bn) 
admits a low energy effective hamiltonian which involves only 
the unlinked atoms (B\, A#), given by 

k* N 




<H eff ~ - 



(t/af 





k N 



(8) 



and hence the new vacancy state found in the ABC trilayer for 
the B1-A3 is directly generalizable to B\-An in the multilayer 
rombohedral compounds. This finding exhausts the possibili- 
ties for the types of vacancy states in multilayer graphenes in 
the two more common stacking. 

A very interesting open issue is the fate of these states for the 
twisted multilayers II 3 2L 1 3 3fl although it is probable that the de- 
scribed features will persist for the vacancies having the same 
neighbours in the twisted superlattice as the ones described 
here. 

The magnetic behaviour associated to the two types of va- 
cancies of the bilayer system was analysed in 112811 . The fully 
localised midgap states arising from the new vacancy states in 
the presence of a gap will give rise to fully localised magnetic 
moments that will play a predominant role in the magnetic be- 
haviour of the samples. 
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